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Recently, a new algorithm for the computation of covariant Lyapunov vectors and of corresponding 
local Lyapunov exponents has become available. Here we study the properties of these still unfamiliar 
quantities for a number of simple models, including an harmonic oscillator coupled to a thermal 
gradient with a two-stage thermostat, which leaves the system ergodic and fully time reversible. 
We explicitly demonstrate how time-reversal invariance affects the perturbation vectors in tangent 
space and the associated local Lyapunov exponents. We also find that the local covariant exponents 
vary discontinuously along directions transverse to the phase flow. 



I. INTRODUCTION 

Recently, many concepts and methods of dynamical 
systems theory have turned out to be very useful for 
the characterization and understanding of physical sys- 
tems in and out of thermodynamic equilibrium. For ex- 
ample, for a class of stationary nonequilibrium systems, 
the spectrum of Lyapunov exponents is a convenient tool 
for studying the collapse of the phase-space probability 
distribution onto fractal measures with an information 
dimension smaller than the dimension of phase space. 
In this case, stationarity is achieved with time-reversible 
thermostats P, [2j. Stationary nonequilibrium systems 
with stochastic thermostats may be formulated along 
similar lines 

The aim of this paper is to apply the hitherto rather 
unfamiliar concept of covariant Lyapunov vectors and 
their associated local Lyapunov exponents to some sim- 
ple and pedagogical systems in equilibrium and in non- 
equilibrium stationary states to sharpen the intuition 
for more demanding applications. The systems studied 
include a harmonic oscillator subjected to a two-stage 
chain of Nose-Hoover-type thermostats with a tempera- 
ture which varies with the position of the particle. 

The paper is organized as follows: In the next sec- 
tion we provide the basic theoretical concepts and def- 
initions required for our numerical work. In particular, 
the covariant vectors and their classical counterparts, the 
Gram-Schmidt vectors, are introduced, and their dynam- 
ical evolution is discussed. Section [TTTJ is devoted to an al- 
ternative differential-equation based method for the evo- 
lution of orthonormal perturbation vectors, which may be 
interpreted as continuous re-orthonormalization. In Sec- 
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tion IIVI we specify the protocol for our numerical work, 
both forward and backward in time. Our main example, 
a doubly-thermostated oscillator in a space-dependent 
thermal field, is treated in various subsections of Section 
fVl Section lVTl is devoted to symplectic systems, with reg- 
ular trajectories on a torus or with chaotic behavior, for 
which the differences of the symmetry properties for the 
local Gram-Schmidt and covariant Lyapunov exponents 
are most pronounced. We conclude in Section I VIII with 
some remarks, which also concern the stationary fluctu- 
ation theorem for thermostated systems. 

II. COVARIANT LYAPUNOV VECTORS AND 
LOCAL LYAPUNOV EXPONENTS 

If T(t) denotes the state of a dynamical system of di- 
mension D, its evolution equations are given by 

r = F(r), (i) 

where F is a (generally nonlinear) vector-valued func- 
tion of dimension D. An arbitrary perturbation vector 
6T(t) in tangent space evolves according to the linearized 
equations 

6T = J(T)6T, (2) 
where the dynamical (Jacobian) matrix J is given by 

OF 

m = w - 

The stability of a trajectory in a D-dimensional phase 
space is determined by a set of D (global) Lyapunov ex- 
ponents, which are the time-averaged logarithmic rates 
of the growth or decay of the norm of some perturba- 
tion vectors, which must be oriented 'properly' in tan- 
gent space at the initial time. Formally, let T(0) denote 
the state of the system at time 0, the state at time t is 
given by T(t) = 0*(r(O)), where the map : T -4 T 



2 



defines the flow in the phase space I\ Similarly, if <5r(0) 
is a vector in the tangent space at the phase point r(0), 
at time t, it becomes ST(t) = D0*|r(o) £r(0), where Dft 
defines the tangent flow. It is represented by a real but 
generally nonsymmetric D x D matrix. The multiplica- 
tive ergodic theorem of Oseledec asserts that there 
exist 'properly oriented' and normalized vectors v e (T(0)) 
in tangent space at t — 0, which evolve according to 



Z^'lrco) v e (r(0))=t/ (T(t)), 



(3) 



and which generate the Lyapunov exponents on the way. 
1 

t^±oo \t 



lim ^ In ||D^| r(0) v e (T(0)) || (4) 



for all £ £ {1, ...,D}, both forward and backward in 
time (for time-reversible systems). (Strictly speaking, 
this formulation is only correct for nondegenerate expo- 
nents Af. If two such exponents become identical, the 
respective vectors must be replaced by a covariant sub- 
space spanned by the vectors. Since in our applications 
below, this happens only for the symplectic systems in 
thermodynamic equilibrium discussed in Sec. IVI1 there 
is no danger of misinterpretation, and we avoid the ad- 
ditional notational complexity. The case of degenerate 
exponents is treated in detail in Ref. Q)- Because of the 
property described by Eq. (|3]), the vectors v l are called 
covariant. Loosely speaking, covariant vectors are co- 
moving (co-rotating in particular) with the tangent flow. 
As will be shown below, this property of co-rotation is 
responsible for the fact that the evolution of their length 
in the forward and backward directions of time (for time- 
reversible systems) is intimately connected, a symmetry 
not enjoyed by other perturbation vectors. For numeri- 
cal reasons, it is still necessary to normalize the vectors 
periodically at times t n = tit, such that Eq. ((4]) becomes 



1 



N-l 

X e = lim — V ln||^ T |r„ v e (T r , 

71 = 



(5) 



where we use the abbreviated notation T(t n ) = T n . 
\\v e (T n )\\ = 1 at the beginning of each interval of length 
r. Generally its norm differs from unity at the end of the 
interval. 

Up to very recently, no practical algorithm for the com- 
putation of the covariant vectors was available. The clas- 
sical algorithm for the computation of Lyapunov expo- 
nents p, E| is based on the fact that almost all volume 
elements of dimension d < D in tangent space (with 
the exception of elements of measure zero) asymptoti- 
cally evolve with an exponential rate, which is equal to 
the sum of the first d Lyapunov exponents. Such a d- 
dimensional subspace may be spanned by d orthonormal 
vectors, which are constructed by the Gram- Schmidt re- 
orthonormalization procedure and, therefore, are referred 
to as Gram-Schmidt (GS) vectors g e (T(t)). The evolu- 
tion during the time interval r = t n — t n —\, 



generates a set of non-orthonormal vectors , 
{g e (r„) , £ — 1, • • • . D}, which after Gram- Schmidt 
re-orthonormalization PJJ, LLl| , 

g 1 (Tn) 

Us 1 (r«) || ' 

9 e (Tn) - Et\ (9 e (Tn) ■ g k (TQ) g k (r n ) 

US' (r n ) - Eti (s e (r„) • g k (r B )) g k (r„) | 



i 1 (r r , 



9 l (r„) 



(where i consecutively assumes the values 1, • • • ,D) be- 
come the orthonormal starting vectors for the next inter- 
val. The vectors g e are not covariant, which means that, 
in general, the vectors are not mapped by the linearized 
dynamics into the GS vectors at the forward images of 
the initial phase-space point [12[ . As a consequence, they 
are also not invariant with respect to the time-reversed 
dynamics. The Lyapunov exponents are computed from 
the normalization factors, 



N-l 



Xf = 



lim 

at->oo Nt 



^£ ln l^( r ")ll> 

n=0 
1 N ^ 

lim V ln||g £ (r n ) 



n=0 



k=l 



{g e (T n )-g k (T n ))g k (T n ) 



(7) 



for £ = 2, • • • , D. 

Recently, a reasonably fast algorithm for the compu- 
tation of covariant Lyapunov vectors was presented by 
Ginelli et al. (l2l |. which first requires the construction 
of the Gram-Schmidt vectors by a forward integration in 
time. In a second step, this stored information is used to 
obtain the covariant vectors by a backward iteration in 
time. For details of this algorithm we refer to their paper 
[lH and to our previous work (tJ. 

A local Lyapunov exponent characterizes the expan- 
sion, or shrinkage, of a particular tangent vector during 
a (short) time interval r. From Eqs. (|5|) and ([7]) local 
exponents for the covariant and Gram-Schmidt vectors 
are obtained for a time t n = nr at the phase point r„: 



A? 0V (O = -hi||£»^|r n _ 1 i/(T n _i)| 
r 



(8) 



for £ = 1, • • • , D, and 



Ap S (t„) 
Ap S (*„) 



in||^(r„) 



J2(9 e (r n )-g k (r n ))g k (T n )\l (9) 



k=l 



Dflr^ 9 e (T n ^)=g e (r. 



(6) 



for I = 2, • • • , D. 
Since the spaces 



(10) 
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are covariant subspaces of the tangent space for all I, we 
have v e (t) € g 1 ^) © • • • © g £ (t). If /3 u (t) denotes the an- 
gle between the respective covariant and Gram- Schmidt 
vectors v it) and g it) at the specified time, the compo- 
nent of the normalized vector v (t n _i) in the direction 
of g e it n -i) is given by cos/3«(i n -i). During r, this vec- 

tor component grows by a factor exp{A^ t}, whereas 
the norm of the vector itself grows by exp{A^ ov T}. At 
the end of the interval, equating the vector components 
of v \tn) in the new direction of the re-orthonormalized 
vector g e it n ), one obtains 



ApS (0 = Ar(fm) + l ln ^fa) 

T COS/%(t„_i) 



(11) 



£ = 1, • • • , £>. This relates the local exponents for the GS 
and covariant vectors. 

If we consider the limit r — > implying continuous 
re-orthonormalization of the and normalization of the 
v , Eq. (jTTJ) becomes 



Ap S (t) = A, cov (t) 



tan/3«(i) 



This is most easily achieved with a matrix of Lagrange 
multipliers constraining the vectors to unit length and 
enforcing orthogonality of the g e [T3T[l5l | . We shall return 
to this point in the following section. 

For time-continuous systems, these relations are gen- 
eral and are not restricted to any particular model. Below 
they will be applied to a variety of models mentioned in 
the introduction. 

The global exponents are the time averages of the local 
exponents over a long trajectory tracing out the whole 
ergodic phase space component, and are the same for the 
covariant and Gram-Schmidt cases, 



N-l 



N-X 



Xe = Km 1 y A™ v it n ) = lim 1 V A 



GS 



n=0 



AT-foo N 



(tn). 



n=0 



Whereas the global Lyapunov exponents do not depend 
on the particular metric and the choice of the coordinate 
system, the local exponents do. This will become appar- 
ent in Section IVII for a scaled harmonic oscillator. For 
particular applications of the local exponents this must 
be kept in mind. 

All systems we consider here are invariant with re- 
spect to time reversal. This property leaves the equa- 
tions of motion in phase and tangent space unchanged if 
the signs of all momentum-like variables and of time are 
reversed, but leaving all position variables unchanged. 
This implies that there exists a smooth isometry / of 
phase space, such that 70* = cf> I. In practice, an in- 
tegration of the equations of motion backward in time is 
carried out with reversed momentum-like variables and 
a positive time step. After reaching the endpoint, the 
signs of all momentum-like variables need to be reversed 
again and the time variable properly adjusted. Alterna- 
tively, and even more easily, the integration of the motion 



equations may proceed without changing the sign of the 
momentum-like variables but with a negative time step. 
There is also no sign change after reaching the end point 
in this case. A comparison of both methods yields iden- 
tical results. Where necessary, we indicate the forward 
and backward directions of time by upper indexes (+) 
and (— ), respectively. If this index is omitted, the for- 
ward direction is implied. 

We have mentioned already that the classical al- 
gorithm invoking Gram-Schmidt re-orthonormalization 
carefully keeps track of the time evolution of d- 
dimensional volume elements, SVd, for any d < D, which 
proceeds according to @, H(| 



dln8V d (t) 
dt 



If the total phase volume is conserved as for symplectic 
systems, the following sum rule for the Gram-Schmidt 
local exponents holds at all times: 



y / A? s (t) = o. 



(12) 



In this symplectic case we can even say more. For each 
positive local GS exponent there is a local negative GS 
exponent such that their pair sum vanishes [171 . 



<+)Ap S (t) 
( ->Ap s (t) 



_(+) A GS 
.(-) A GS 

1Y D+X- 



>(t), 
!(*)• 



(13) 
(14) 



As indicated, such a symplectic local pairing symmetry 
exists both forward and backward in time. But, gener- 
ally, the GS local exponents do not show the symmetry 
with respect to time-reversal invariance. Thus, 



{ - ] A? S it) ^ - (+) A 



GS 

D+X- 



(15) 



No such symmetries exist for non-symplectic systems. 
Examples are provided below. 

The situation is very different for the covariant local 
Lyapunov exponents. In their case, the vectors are still 
re-normalized, but the angles between them remain un- 
changed, which effectively destroys all information con- 
cerning the d-dimensional volume elements. Thus, no 
symmetries analogous to Eqs. ([13]) and ([T4|) exist. In- 
stead, the re-normalized covariant vectors faithfully pre- 
serve the time-reversal symmetry of the equations of mo- 
tion, which is reflected by 

H ArW = - (+) A C „° + V i-, it) for I = 1, ... , D, (16) 

regardless, whether the system is symplectic or not. This 
means that an expanding co-moving direction is con- 
verted into a collapsing co-moving direction by an ap- 
plication of the time-reversal operation. Of course, the 
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forward and backward local exponents in Eq. II 6|) refer 
to the same phase space point T(t). 

These symmetry properties may be considered the 
main conceptual differences between the Gram-Schmidt 
and covariant viewpoints. 

Before leaving this section, a short remark concern- 
ing the commonly-used term "time-dependent exponent" 
seems in order. Primarily, this quantity is a function of 
the phase point and should only be called a "local" expo- 
nent. Its value may be different whether the phase point 
is reached from the past, forward in time (+), or from 
the future, backward in time (— ). 



III. DIFFERENTIAL APPROACH TO LOCAL 
LYAPUNOV EXPONENTS 

Equation @ precisely reflects the numerical procedure 
for the computation of local GS exponents for finite time 
intervals r. But it is also possible to obtain a differential 
version for r — >• 0. Goldhirsch et al. derived a full set 
of differential equations for the Gram-Schmidt vectors g l 
131, 



g 1 = Jg 1 -R u g 1 , 



(17) 



g e = Jg e - Ru 9 e -Y, V**™ + 9 m , (18) 



where in the last equation I = 2, • ■ • , D. We have demon- 
strated [1J, ll5| that the matrix elements 



Rim (r(t)) = (g'fjg 



(19) 



may be understood as Lagrange multipliers enforcing the 
orthonormalization constraints g e ■ g m — 5g m (equal to 
unity for £ = m, and zero otherwise). Here T means 
transposition. Most importantly, the diagonal elements 
are the local Gram-Schmidt exponents: 



Ap s (r(t)) = R u (T(t)) = (g e ) T Jg e 



(20) 



This expression nicely underlines the local nature of the 
exponents. 

We have verified for the doubly-thermostated heat con- 
duction model discussed in Sec. IVl below that the direct 
integration of the Eqs. (|17I18|) provide local GS expo- 
nents according to Eq. (|20p . which agree extremely well 
with the results obtained from a direct application of the 
GS algorithm, Eq. ©, for a reasonably-small time inter- 
val r. This agreement also persists for the time-reversed 
dynamics. 



IV. NUMERICAL CONSIDERATIONS 

In this section we remark on a few aspects of our im- 
plementation of the algorithm for the computation of the 



covariant Lyapunov vectors, which we apply in the fol- 
lowing sections. Reduced units are used for the vari- 
ous models treated below. For convenience, we specify 
already here the adopted values (in reduced units) for 
some time parameters: t u = 6 x I0 4 , t a — 5 x 10 4 , 
and to — 100. Their meaning is explained below. For 
the integration of the equations of motion, a 4th-order 
Runge-Kutta algorithm with a time step dt — 0.001 is 
used. This time step is chosen such that the trajectory 
is correct to double-precision accuracy. For the interval 
between successive Gram-Schmidt re-orthonormalization 
steps - respective covariant vector normalizations - we 
choose t — lOdi = 0.01. This number is a (very conser- 
vative) compromise between the achieved reduction in 
storage requirements as outlined below, and the preci- 
sion of integration forward and backward over the same 
interval. The time to is chosen such that in the inter- 
val —to < t < to accurate Gram-Schmidt and covariant 
Lyapunov vectors are available. 

The simulations are carried out with the following pro- 
tocol: 

Phase 1 (forward integration from — t u to +t ul ): 
Starting with arbitrary initial conditions at a time — t u 
and using a positive integration time step dt > 0, the evo- 
lution of the reference trajectory T(t) and of the full set of 
Gram-Schmidt vectors is computed in the forward direc- 
tion of time up to a time +t u . The reference trajectory 
and the Gram-Schmidt vectors are stored for every 10 
time steps, 10e?i = r, along the way. The Gram-Schmidt 
vectors are used in phase 2 for the construction of the co- 
variant vectors, and the reference trajectory is required in 
phase 3 for the computation of the time-reversed Gram- 
Schmidt vectors. The Lyapunov spectrum {( + )A^} is 
accumulated for times —t a <t<t u , for which the ori- 
entations of the Gram-Schmidt vectors are fully relaxed. 
Phase 2 (backward iteration from t u to — to): Start- 
ing at t u , the covariant vectors are computed by iter- 
ating back to a time —to- The details of this algo- 
rithm are given in Ref. [7]. Since the forward GS- 
vectors, stored during phase 1, are now used in reversed 
order, the consecutive order of the covariant vectors 
• • • , v (t n ), v e (t n -i), ■ ■ ■ has to be reversed for the com- 
putation of the corresponding local exponent, 

(+) A COV r , , 1 \\vHtn)\\ 

{tn) ~ r \y(t n ^w 

or, alternatively, the sign of the local exponents must 
be reversed. The time averaging for the global Lyapunov 
spectrum {^A£° v } is carried out for times t a > t > —to. 

The following two phases are only required for the com- 
putation of the local time-reversed Gram-Schmidt and 
covariant exponents. 

Phase 3 (backward integration from t u to —tu): 
With arbitrary initial conditions at time t u , the Gram- 
Schmidt tangent space dynamics is followed backward in 
time up to —t u . To counteract the Lyapunov instabil- 
ity, it is essential for this computation to use the same 
reference trajectory stored in phase 1, where the sign of 
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the momentum-like variables ( p, z, and x for the doubly- 
thermostated oscillator) is left unchanged, but with the 
time step reversed to —0.001. For an accurate computa- 
tion of the backward GS vectors, the reference trajectory 
at every integration step is required. Since in phase 1 
this information was stored for only every 10th step (to 
save computer storage) , the forward reference trajectory 
is piece- wise re-computed with stored phase-space points 
as initial conditions. The backward Gram-Schmidt vec- 
tors are again stored for every 10th time step replacing 
the forward vectors of phase 1. If time is reversed, the 
stable directions become unstable and vice versa. The 
Lyapunov spectrum {("•'Ap^} is accumulated in the in- 
terval t a > t > —t u . 

Phase 4 (forward iteration from — t u to +to): Anal- 
ogous to phase 2, in this final stage the covariant Lya- 
punov vectors for the reversed time direction are com- 
puted with the help of the time-reversed Gram- Schmidt 
vectors from phase 3. The respective Lyapunov spectrum 
'A^ ov } is accumulated for times — t a < t < to- 
It may be noticed that in the interval —to<t< +tg all 
local properties are available with the Gram-Schmidt and 
covariant vectors fully relaxed both forward and back- 
ward in time. Therefore, the detailed analysis of local 
(time-dependent) Lyapunov exponents in the following 
sections is carried out in this regime. 



V. DOUBLY-THERMOSTATED OSCILLATOR 
A. Description of the model 

Here we consider a simple model which already has 
many ingredients in common with much more involved 
physical systems. It exhibits chaotic equilibrium and 
stationary nonequilibrium states and collapses onto a 
limit cycle for very strong driving. It consists of a one- 
dimensional harmonic oscillator, which is coupled to two 
consecutive stages of Nose-Hoover thermostats with a 
space-dependent temperature T(q). The equilibrium ver- 
sion of this model was first considered by Martyna, Klein 
and Tuckerman [l8j . Its nonequilibrium properties were 
consecutively studied by us in some detail [l^, HO] , but 
without considering covariant vectors. This paper is also 
intended to augment this work correspondingly. 

The equations of motion, expanded with two thermo- 
stat variables z and x, are (19L l20l| 

4 = P, 

p = -q- zp, 

z = p 2 — T(q) — zx, 

x = z 2 -T(q), 

where the position dependent temperature is given by 
T(q) = l + £tanh(g). 



The control parameter e coincides with the temperature 
gradient at q = 0. These equations are written in the 
most simple reduced form with all arbitrary parameters 
of the model set equal to unity. The system is not sym- 
plectic. On average, the oscillator picks up energy from 
the thermostat whenever it is in a region of high temper- 
ature (q > 0), and releases it again in low-temperature 
regions (q < 0). 

B. Global properties 

For a typical non-equilibrium state, (e = 0.25), the 
global Lyapunov spectrum was computed by four inde- 
pendent methods, applying the protocol outlined in Sec. 

El 

Phase 1 : GS exponents in forward direction of time, 
|(+) A GS} = {o.053i, 0.0000i, -0.034 4 , -O.O867}, 

Phase 2 : covariant exponents in forward direction of 
time 

|(+) A cov| = {0.053 6 ,0.0000i,-0.035i,-0.086 2 }, 

Phase 3 : GS exponents in backward direction of time 
{(-)Ap S } = {0.086 7 ,0.034 4 ,0.0000 3 ,-0.053i}, 

Phase 4 : covariant exponents in backward direction of 
time 

|(-) A cov| = {0.087i,0.033 7 ,0.00000i,-0.052 5 }. 

The last digit of each number is rounded accordingly. 
Considering the smallness of the exponents and the 
rather involved numerical procedures, the agreement be- 
tween the independently-determined global spectra is 
very satisfactory. A comparison of the forward and back- 
ward dynamics reveals the theoretically expected symme- 
try for the global Lyapunov exponents [!, [2l[ , 

1-1 \t = -^Xo+i-i for I = 1, ■ • • ,D. (21) 

If the temperature gradient e is varied over a wide 
range, significant changes of the spectrum become evi- 
dent. This is shown in the top panel of Fig. [TJ There 
exist a number of distinct regimes with different qualita- 
tive behavior. 

For e < 0.18, the spectrum changes but little with 
e, and the Kaplan- Yorke dimension is only weakly re- 
duced with respect to the full phase space dimension, as 
is shown in the bottom panel of Fig. [TJ The dissipation 
due to the weak heat current influences the appearance 
of the chaotic phase-space trajectory very little. An ex- 
ample of a projection of such a trajectory onto the qpz- 
subspace is provided in the top panel of Fig. [5] 

For 0.18 < £ < 0.26, the trajectory is more and more 
attracted to a weakly unstable periodic orbit (see the 
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FIG. 1: (Color online) Temperature-gradient dependence of 
all four Lyapunov exponents (top panel) and of the Kaplan- 
Yorke dimension (bottom panel) for the doubly-thermostated 
oscillator. For e > e c = 0.26314, the trajectory collapses onto 
a limit cycle. 
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bottom panel of Fig. [5]), which for e c 0.26312 turns 
into a stable limit cycle as shown in the top panel of 
Fig. [4] The nature of this transition may be established 
by considering the Floquet multipliers \xi, t = 1, • • ■ ,4 
for the fixed points of the Poincare map, defined by q = 
0, for e > e c ~ 0.26312. Whereas fix = 1 and /i4 < 
0, a single mutiplier /i2 crosses the unit circle on the 
real axis at the point A corresponding to e c in Figure 
[3J Such a behavior is characteristic of a period doubling 
bifurcation [22j], where, possibly, the chaotic attractor 
disappears in a boundary crisis bifurcation. This point 
will be studied separately (23|. Increasing e further, the 
Floquet multipliers )i2,3 vary as indicated by the arrows 
and become complex conjugate to each other for e w 
0.26319 (point B in Fig. 



For e w 0.417, there is another transition changing 
the two-loop limit cycle into a single-loop orbit. This 
is illustrated in the bottom panel of Fig. @] and will be 
studied separately 23 [. 
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FIG. 2: (Color online) Floquet multipliers fi2 and [13 in the 
complex plain for e > 0.26312. UC denotes the unit circle, 
and A is the bifurcation point. The multipliers are complex 
conjugate to each other for e > 0.26319 as indicated by B. 



C. Local Lyapunov exponents 

In Fig. [S] we apply Eq. (fTT| to the doubly-thermo- 
stated oscillator in a stationary chaotic state, e = 0.25. 
For I = 1 the respective time-dependent exponents are 
identical and are not shown. The case i = 2 is treated in 
the figure. The dashed green line denotes the covariant 
local exponent, the smooth red line is for the local GS- 
cxponents, which is directly obtained from the simula- 
tion invoking Gram-Schmidt re-orthonormalization. The 

time interval r is 0.01. The blue points for Ap S (i) are 
computed with Eq. (|11[) . where the covariant exponent 
A^ ov (t) and the angle Pu(t) are taken from the simula- 
tion. The agreement is convincing. Similar results are 
also obtained for t — 3 and 4 (not shown). 

In the bottom panel of Fig. [S]we demonstrate, for I — 
2, the general time-reversal symmetry for the local (time 
dependent) covariant exponents (see Eq. (fTTJ)) which also 
gives rise to the symmetry of the global (time-averaged) 
exponents already encountered in Eq. (f2"Tj) . For 1=1 
the symmetry is also fully obeyed but not shown. 

As emphasized already in Eq. (fl5)) . the local Gram- 
Schmidt exponents generally do not have this symme- 
try. This is explicitly shown in the top panel of Fig. 
[SI See also Ref. (2(j, where the same observation was 
made. Only the subspaces in Eq. (fT0|) spanned by con- 
secutive Gram-Schmidt vectors have a simple dynamical 
interpretation, but not the GS-vectors themselves. The 
orthonormal GS-vectors are oriented such that for the 
tangent space, tangent to the phase flow at the phase 
point T(i), the subspaces (t) 8 ••• © ( ~ ) g e , with 

I E {1, ■ • • , D}, are the most unstable subspaces of di- 
mension I going from time t to —00 (i.e. the most stable 
subspaces of dimension I in the future). Although time 
reversal converts a most stable subspace of dimension i 
into the most unstable subspace with the same dimen- 
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FIG. 3: (Color online) Projection of a short chaotic trajec- 
tory for e — 0.10 (top) and e = 0.25 (bottom) onto the qpz- 
subspace. 

sion, and vice versa, there is no obvious correlation of 
the instantaneous Lyapunov exponents ' - )Ap^(t) and 
WAg^) brt = l,...,D. 

It is interesting to follow the time dependence of the 
covariant local exponents, or more correctly expressed, 
their variation for consecutive state points along the 
phase space trajectory (see Fig. [7]). One observes that 
the order of the exponents fluctuates and may even be 
totally reversed with Al° v (t) being most negative and 
A4° v (i) most positive. Also the dimension of the sta- 
ble and unstable manifolds changes along the trajectory. 
This indicates that the system is far from being hyper- 
bolic. We address this point more closely in the following 
subsection. 



D. Hyperbolicity 

We infer from Eq. (TTTI) that the difference between 
the local covariant and Gram-Schmidt exponents stems 
from the fact that the angle between the respective vec- 
tors deviates significantly from zero and varies with time. 
But also the angles between covariant vectors, ctij(t) = 
arccos[(v l • ■w-')/|t7 l ||t7 3 '|] significantly change with time. 



e = 0.28 




FIG. 4: (Color online) Projection of the limit cycle for e = 
0.28 (top) and e = 0.45 (bottom) onto the gpz-subspace. 



This is demonstrated in the bottom panel of Fig. [5] for 
the same nonequilibrium state (e = 0.25) of the doubly- 
thermostated oscillator discussed previously. There is 
an intermittent tendency of any two pairs of vectors to 
get parallel or antiparallel to each other. The proba- 
bility distributions for these angles, 7t(ay), are shown 
in the top panel of Fig. [5] and confirm this observa- 
tion. Although the angles a do not become strictly 
zero - the vectors could not separate anymore after such 
an event, which is not observed - the large probabil- 
ity for angles close to zero or ir is noticeable. As was 
mentioned before, the associated local covariant expo- 
nents are out of order for most of the time as in Fig. 
[7] If Pi denotes the probability for Af ov to be out of 
order with respect to any of the other exponents, one 
finds for the doubly-thermostated oscillator (e = 0.25) 
{P u ... ,P 4 } = {0.650,0.813,0.840,0.645}. This clearly 
demonstrates the strong entanglement between the co- 
variant vectors. If the local exponents are time averaged 
along the trajectory for time intervals A, the analogous 

probabilities pf for the time-averaged exponents A^ ov 

scale according to P cx A _7i with 7 > for large- 
enough A. This shows that the domination of the Os- 
eledec splitting is violated for finite times. 

Such a behavior is in contrast to the covariant dynam- 
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FIG. 5: (Color online) Time-dependence of the local Lya- 
punov exponents A2 (t) in the forward direction of time for the 
doubly-thermostated oscillator with e = 0.25. The smooth 
red curve denotes Gram-Schmidt exponents, directly obtained 
with a re-orthonormalization procedure, the blue points are 
computed with Eq. (jlll) . using the covariant time-dependent 
exponents ov (also shown as green dashed line) and the 
angles /3u as input. For clarity, only every 20th point is de- 
picted. 



ics of hard-disk systems, for which the covariant vectors 
tend to avoid becoming parallel or antiparallel Q- Thus, 
whereas the hard-disk system is hyperbolic, the doubly- 
thermostated oscillator is not. 



E. Singularities of the local Lyapunov exponents 

In the direction of the flow, the local Lyapunov ex- 
ponents clearly are smooth functions of the time and, 
hence, of the phase-space position along the trajectory, 
see Fig. [5] But transverse to the flow this need not be the 
case. Indeed, for the periodic Lorentz gas it was noted 
by Gaspard [H], HE] that the local stretching factors are 
discontinuous transverse to the flow. Since this model in- 
volves hard elastic collisions of point particles with space- 
fixed scatterers, the observed discontinuity might still be 
thought to be a consequence of the discontinuous na- 
ture of the flow. However, Dellago and Hoover showed 
26] that this is not the case. They found a discontinu- 
ous local exponent Af" 3 along a path transverse to the 
flow even for a time-continuous Hamiltonian system, a 
chaotic pendulum on a spring. Of course, their result 
also applies to A^ ov for that model. Here we provide 
evidence for the doubly-thermostated oscillator in equi- 
librium (e = 0) that all local covariant exponents are 
discontinuous along directions transverse to the flow. 

For this simulation we slightly modify the protocol of 
Section El 

Phase 0: Starting at a phase point T s at time zero, 
the reference trajectory is followed backward in time to 
— t u — —60,000 and is periodically stored for intervals r 
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FIG. 6: (Color online) Doubly-thermostated oscillator for 
e = 0.25. Top panel: The Gram-Schmidt local Lyapunov 
exponents do not display time-reversal symmetry. Bottom 
panel: Display of time-reversal symmetry by the covariant lo- 
cal exponents, (+) Af ov = - (-) Ag^-!! for e = 2 > Analogous 
curves are obtained for the other £, but are not shown. 



along the way. 

Phase 1: The next phase is identical to phase 1 of Sec- 
tion |IV] with one essential difference: For — t u < t < 0, 
the previously-stored reference trajectory is now used in 
the forward direction of time for the computation of the 
Gram-Schmidt vectors, which assures that the trajectory 
precisely arrives at r s at time zero in spite of the inher- 
ent Lyapunov instability. For < t < t u the simulation 
proceeds as in phase 1 of Section HVl 
Phase 2: This is identical to phase 2 of Section HVl and 
provides us with the covariant vectors and the respective 
local exponents in the interval —to < t < to and at the 
time t = in particular, when the state coincides with 
the selected phase point r s . 

The whole procedure is repeated for starting points 
r s = To + s x (0, 0, 1, 0) on a straight line parallel to the z- 
axis, which is parametrized by s. This line is transversal 
to the flow, as may be expected from Fig. [5] 

As an example, we plot in the top panel of Fig. [5] the 
local covariant exponent A4° v (i) as a function of time 
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FIG. 7: (Color online) Doubly-thermostated oscillator with 
e = 0.25: Time dependence of all four local covariant Lya- 
punov exponents. 



for 500 initial points T s separated by As = 1 x 10~ 4 . It 
should be noted that the scale on the t axis, converted 
into distances along the trajectory, is about 200 times 
coarser than that on the s axis. The red line for t = 
connects the local exponents for the selected initial states 
r s . This curve for A2° v (t) is also reproduced in the bot- 
tom panel of Fig. |H] together with an analogous result 
for Ai° v (t). Both curves exhibit singularities on many 
scales showing singular fractal character. There are no 
obvious correlations between the two curves. The singu- 
larities are due to bifurcations in the past history of the 
trajectory. In view of Fig. [SJ such a bifurcation may be 
visualized, for example, by a transition of the trajectory 
from the neighborhood of an unstable periodic orbit to 
the neighborhood of another with a different number of 
loops. 

One may raise the question (as has been done by one 
of the referees) , how reliable the curves in Fig. [9] are in 
view of the chaotic nature of the flow and problems of 
shadowing due to the finite computational accuracy. An 
increase of the Runge-Kutta integration time step dt by 
a factor of four has no noticeable effect (less than 0.1%) 
in Fig. [§J which also proved completely insensitive to a 
reduction of the relaxation time t u of the algorithm by 
a factor of two and of an increase of the time r between 
successive re-orthonormalization steps by the same fac- 
tor. This robustness, however, does not apply to the lo- 
cal exponents A2° v and A3° v (not shown), which belong 
to the two-dimensional central manifold for this equilib- 
rium system. The respective covariant vectors span this 
subspace, but their precise orientations and their local 
exponents are affected by details of the algorithm and do 
not have direct physical significance. 

For nonequilibrium stationary states the singular char- 
acter of the local exponents in transverse directions is 
expected to be even more pronounced, since even the 
phase-spa ce p robability density becomes a multifractal 
object [l|, [l4(. For the covariant exponent this cannot 
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FIG. 8: (Color online) Doubly-thermostated oscillator with 
e = 0.25. Top panel: Probability distributions for the an- 
gles ctij between covariant vector pairs specified by the labels. 
Bottom panel: Time evolution of the angles a%j between the 
covariant vectors v 1 and v 3 . 



be shown with the present algorithm. The reason is 
that during the time-reversed simulation in phase 0, the 
phase volumes collapse yielding negative Lyapunov ex- 
ponent sums. Since in phase 1 this trajectory is followed 
in the opposite direction, the respective phase volumes 
expand providing a positive sum of Lyapunov exponents, 
but only up to time zero. For positive times the refer- 
ence trajectory is calculated anew from the motion equa- 
tions, again yielding contracting phase volumes. Thus, 
the character of the flow changes at t = and the Gram- 
Schmidt vectors at first are non-relaxed and point into 
wrong directions for positive times. Since these vectors 
are required for the computation of covariant vectors at 
and near zero time, the algorithm cannot be used to ob- 
tain the covariant vectors and respective local exponents 
at a pre-determined point T s in phase space. For equilib- 
rium states this restriction does not apply and the local 
exponents may be computed for pre-specified phase-space 
points. 
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FIG. 9: (Color online) Covariant local Lyapunov exponent, 
A4° v , for the doubly-thermostated oscillator in equilibrium. 
Top panel: Time dependence along trajectories, which, for 
t = 0, are at the specified initial points T s introduced in the 
main text. These phase points lie on a straight line transverse 
to the flow, and s specifies the precise position. The variation 
of the local exponent along this straight line is shown as a red 
curve. Bottom panel: The red curve is a magnified view of 
the red line of the upper panel, representing the variation of 
A4° v along a straight line transverse to the flow in the phase 
space. The blue line is an analogous curve for Ai° v . 



VI. LOCAL LYAPUNOV EXPONENTS FOR 
SYMPLECTIC SYSTEMS 

So far we have omitted to mention that we use a par- 
ticular metric in phase space. Whereas the global expo- 
nents are independent of the metric, the local exponents, 
covariant or Gram-Schmidt, clearly are not. We demon- 
strate this with the most sim ple symplectic example, a 
scaled harmonic oscillator 1 23, u3 with Hamiltonian 



H s {p,q) = - 



and equations of motion 



where s is a scaling parameter. For the 'natural' scaling, 
s = 1, the global and local exponents vanish. But for the 
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FIG. 10: (Color online) Scaled harmonic oscillator with s = 2: 
Time dependence of the local Lyapunov exponents. Top 
panel: The smooth red lines are local GS exponents for I — 1 
and 2 as indicated by the labels. The results from direct sim- 
ulation and from theory are undistinguishable on this scale. 
The points indicate the reconstruction of with the 

help of Eq. (|lip . where ' + ^A2° V is shown in the lower panel. 
Note that the GS and covariant exponents for £ — 1 are iden- 
tical. Bottom panel: Demonstration of the time-reversal in- 
variance for the covariant exponents. 

scaled case, s ^ 1, the local Lyapunov exponents do not. 
They depend on the metric and, for that matter, on the 
choice of the coordinate system, be it Cartesian or polar. 

Let us look at this model in a little more detail, since 
the dynamical matrix 



J 











q = s p , p 



-s 2 q, 



(22) 



for this linear model does not depend on the phase point 
and allows for a complete analytical solution [28( for the 
tangent vector dynamics. Still, the model is a bit pe- 
culiar, since there is no global exponential ordering of 
tangent vectors familiar from the Gram-Schmidt algo- 
rithm, and the considerations of Sec. II VI lose their mean- 
ing. Any (unit) vector, with arbitrary initial condition 
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(phase), which is a solution of Eq. (|17j) . may be taken 
as the first Gram-Schmidt vector (or covariant vector for 
that matter). We fix this arbitrary phase by requiring 
that g 1 coincides with the normalized phase-space veloc- 
ity, which is associated with a vanishing exponent and is 
a solution of Eq. (|17p . as may be explicitly shown. Thus, 



8qi 
Spi 



-1/2 



p/s 2 
— s 2 q 



(23) 



where we denote the perturbation components of g 1 by 
Sqi and Spi). From the constrained motion equation (|17[) 
one obtains 



Sq 1 = s 2 6pi - RnSqx 
Sp 1 = -s 2 5qi - RnSpi. 



(24) 
(25) 



Since g 1 is constrained to unit length, Sq\6q 1 + 5p\5p 1 = 
0, and noting that Aj-" 3 = Ru, we obtain from these 
equations 



GS 



SqiSpi 



(26) 



Insertion of 5q\ and Spi from Eq. (|23p yields an analytic 
expression for the local Gram-Schmidt exponent A^ as 
a function of the phase-space point (q,p). The second 
Gram-Schmidt vector is perpendicular to the first, and 
its associated GS exponents immediately follow from the 
conservation of phase-space volume: 



Ap S (g,p) = -A?°(q >P ). 

In the top panel of Fig. [TU] the local Gram-Schmidt 
exponents for the scaled harmonic oscillator for s = 2 
are shown as a function of time. The initial conditions 
q(0) = and p(0) = 1 are used, for which the solution of 
Eq. (f2"2"|) becomes 



GS/ 



9(t) 



' sint , p{t) = cost. 



Both computer simulation results and the theoretical ex- 
pressions for A^ and A^ are shown, which agree so 
well that they cannot be distinguished and appear only 
as a single smooth red line in the figure. 

The upper and lower bounds of the local exponents are 
also easily obtained from Eq. (|2"6"|) 28] , 



A 



mm, max 



+ 2 



2 )/2- 



v 1 and v 2 . In the lower panel of Fig. [T0]the covariant ex- 
ponent A2° v (i) is shown as it is obtained from the simu- 
lation. If this function is used to reconstruct Ap^(<) w ith 
the help of Eq. (1X11) , the (blue) dots in the upper panel 
of Fig. [TUlare obtained, where cos(/?22) = g 2 • v 2 is also 
taken from the simulation. The agreement is very good. 
In the lower panel of Fig. [TU]we also plot ^A^ m {t) for 
the time-reversed dynamics. The time-reversal symme- 
try of Eq. (HD, 

(-)A? ov (i) = -(+)AC° V (<), 

is nicely displayed. 

As a slightly more involved example, we compute the 
four local Lyapunov exponents for the symplectic Henon- 
Heiles system [29j with Hamiltonian 



H 



1 



(pl+pl) + l{x 2 + y 2 )+x 2 y 



1 



(27) 



For an energy H — 1 /6, the system is known to be chaotic 
(with a Lyapunov spectrum {0.1277,0,0,-0.1277}), 
where the trajectory visits most of the accessible phase 
space [13,|3l|. Using the protocol of Section [TV] we com- 
pute the GS and covariant exponents and present some 
of the results in Fig. QTJ In the top panel the symplec- 
tic local pairing symmetry of Eq. (|T3")) for (+'Ap^ and 
( + ) A GS 

is shown by the smooth red lines (similar to the 
results of Ref. 3li])- The green dashed line refers to the 

time-reversed exponent (~)Ap^ and clearly emphasizes 
the lack of any time-reversal symmetry as expressed by 
the inequality flTSj) . On the other hand, for the covariant 
exponents precisely this symmetry is evident from the 
lower panel of Fig. [TT1 

To avoid confusion, we note that the 'detailed balance 
symmetry' introduced in Ref. [3l| is not connected with 
the symplectic local pairing symmetry considered here. 
The former only means that for a global exponent to 
become zero, the positive and negative parts of the re- 
spective local exponent along the trajectory must cancel 
each other when integrated over time. 



The extrema are attained whenever the components 5q 
and Sp of the respective perturbation vectors contribute 
equally, Sp — ±Sq. For s — 2 we find A m j n max = 
±1.875 in full agreement with Fig. [TOJ 

The undetermined phases we encountered with the 
Gram-Schmidt vectors also carry over to the covariant 
vectors. But since the latter are computed with the help 
of the former, the choice of phase for g 1 also fixes that for 
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FIG. 11: (Color online) Local Lyapunov exponents for the 
Henon-Heiles system with energy H = 1/6. Top panel: Illus- 
tration of the symplectic local pairing symmetry, Eq. (|13|) . for 

the Gram-Schmidt exponents (+) Ap^ and (+) Ap° (smooth 
red lines). Here, D = 4. Also the inequality Eq. (|15p applies 

i— \ CS 

(£ — 2), as the dashed green line for 1 'A3 certifies. Bottom 
panel: Verification of the time-reversal invariance property 
(|16|) for the covariant vectors specified. 



VII. CONCLUDING REMARKS 

For the doubly-thermostated oscillator in a nonequilib- 
rium stationary state, there is a single vanishing global 
exponent, A2, due to the time-translation invariance of 
the equations of motion. The corresponding covariant 
vector, v 2 (t), needs to be parallel (or antiparallel) to the 
phase-space velocity T(t) = {q(t),p(t), z(t),x(t)}. We 
have verified in our simulation that this is indeed the 
case. The remaining vectors v l , v 3 and v 4 are oriented 
with angles fluctuating between and n with respect to 
r(t). The Gram-Schmidt vectors behave very differently. 
Whereas the vector g 1 is identical to v 1 , the vector g 2 
is not parallel to T(t). Instead, the vectors g 3 and g 4 



are perpendicular to T(t) as expected in view of the co- 
variant subspaces of Eq. (fTU|) . These observations serve 
as convenient consistency checks for the numerical pro- 
cedure. 

One of the remarkable features of the covariant local 
Lyapunov exponents A cov (r(i)) is their singular behav- 
ior transverse to the phase flow, whereas they are ab- 
solutely continuous in the direction of the flow. Fig. [5] 
provides an illuminating example. The singularities are 
consequences of bifurcations in the past history. Still, the 
local exponents are point functions in the phase space 
in the sense that one always gets the same value at the 
state point in question, as long as the trajectory has been 
followed from far enough back and has experienced the 
same history. Due to the uniqueness of the solutions of 
differential equations there is only this path to the state 
point in question. The global exponents, however, are 
time averages of the local exponents along an (ergodic) 
trajectory. 

A final remark concerns the doubly-thermostated 
driven oscillator again. In a driven system (in our case a 
single particle in a non-homogenous thermal field) heat 
and, hence, entropy is generated, which needs to be com- 
pensated by a negative entropy production in the ther- 
mostat to achieve a stationary state. The excess heat is 
transferred from the system to the thermostat (by the 
positive friction zp > 0), where it disappears. It follows 
from the thermostated motion equations in Sec. IV Al that 
the external entropy production (of the reservoir) is given 
by 

d ■ 

S/k = — - T — z + x, 

where k is the Boltzmann constant. In the non- 
equilibrium situation, a full time average (z + x) is nec- 
essarily positive. However, we have verified by simula- 
tion that finite time averages of this quantity numerically 
obey the steady-state fluctuation theorem originally dis- 
covered by Evans, Cohen and Morriss |32| ]. This theorem 
was given a firm theoretical basis by Gallavotti and Co- 
hen [H, H3|, by invoking the so-called 'chaotic hypothe- 
sis' for Anosov-like systems. Although our system is not 
Anosov-like, it still obeys the theorem. 
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